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Molecular-dynamics simulations of KNbOs reveal preformed dynamic chain-like structures, present 
even in the paraelectric phase, that are related to the softening of phonon branches over large regions 
of the Brillouin zone. The phase sequence of ferroelectric transitions is correctly reproduced, showing 
that the first-principles effective Hamiltonian used in the simulations captures the essential behavior 
of the microscopic fluctuations driving the transitions. Real space chains provide a framework for 
understanding both the displacive and order-disorder characteristics of these phase transitions. 
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There are many experimental indications in perovskite 
ferroelectrics that the actual atomic structure in some of 
the phases may be significantly different locally than is 
indicated by the average crystallographic structure de- 
duced from elastic X-ray and neutron scattering. Among 
the earliest of these, diffuse X-ray scattering measure- 
ments of Comes et at, |]|| on KNb0 3 and BaTi0 3 
revealed temperature-dependent streak patterns, which 
were interpreted as evidence for randomly ordered static 
chains (oriented along [100] directions) of distorted prim- 
itive cells, even in the paraelectric cubic phase. Other 
observations such as quasi-elastic central peaks in neu- 
tron scattering H and Raman spectroscopy ||] above 
the phase transition temperature are also indicative of 
preformed clusters of the low temperature phase. More 
recently, pair distribution functions obtained from neu- 
tron scattering measurements up to very high momen- 
tum transfers Q and XAFS measurements || indicate 
the presence of short-range order. First-principles calcu- 
lations provide a powerful tool for probing the structural 
energetics of local distortions. For KNbOs, LAPW lin- 
ear response calculation of the zero-temperature phonon 
dispersion in the cubic perovskite structure reveals its 
instability against the formation of localized chain dis- 
tortions. @ In this paper, we explore the implications of 
these results for dynamical behavior at nonzero tempera- 
ture through the molecular dynamics (MD) simulation of 
an effective Hamiltonian constructed from first-principles 
calculations, establishing the presence and dynamic na- 
ture of localized chain distortions above T c . 

While first-principles calculations have provided a 
great deal of information, they are limited to zero tem- 
perature and to relatively small simulation cells. The use 
of effective Hamiltonians, H e g , to extend the reach of the 



first-principles results has proved to be a useful strategy 
for the quantitative analysis of temperature-driven struc- 
tural phase transitions in real materials. |8|||[l0| These 
Hamiltonians act in the subspace of the full ionic con- 
figuration space which contains the degrees of freedom 
relevant to the transition(s). These include, in partic- 
ular, the "soft mode," identified as the unstable mode 
of the high-symmetry structure which freezes in to pro- 
duce the low-symmetry phase(s). The coefficients in a 
Taylor expansion around the high-symmetry structure 
of the Born-Oppenheimer surface in this subspace are 
determined from first-principles total-energy and linear- 
response results. Nonzero temperature simulations using 
H e g quantitatively reproduce the structural transition 
properties of the full system. 

We constructed H e g for KNbC-3 using the lattice- 
Wannier-function (LWF) method. |J Full details of the 
construction are presented elsewhere, [ fj"l"[ and we give 
only a brief description here. The effective Hamilto- 
nian subspace is defined using a basis of localized and 
symmetrized atomic displacement patterns, called lattice 
Wannier functions, which are constructed to reproduce 
the first-principles unstable polar Ti5 phonon as well as 
unstable transverse optic phonon eigenvectors and fre- 
quencies at other high-symmetry points in the BZ [Q. 
For KNbC-3, the subspace is spanned by one vector de- 
gree of freedom per unit cell, £i a , representing the LWF 
coordinates, where i = unit cell index and a — x,y,z. 
We include as additional degrees of freedom the homoge- 
neous strain tensor, which describes changes in the over- 
all volume and shape of the simulation cell. 

In the LWF basis, the kinetic energy retains a simple 
diagonal form. The potential energy is expressed as a 
Taylor expansion in the LWF coordinates £j a and can be 
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organized as follows: 

U Ij T on — site U ^ short — range Unipolar 

+ Ulwf- strain lastic ■ (i) 

We include anharmonic terms only in the on-site in- 
teraction Uon-site = <f + 5£f + l{£,f x (i y + c.p.), and 
in the lowest order coupling between LWF coordi- 
nates (on-site quadratic) and homogeneous strain (lin- 
ear) V lwf -strain- Strain coupling of this type has been 
shown to be crucial in obtaining the experimentally ob- 
served sequence of ferroelectric phase transitions in per- 
ovskite ferroelectrics. ||[l2| The coefficients of the on- 
site anharmonic and strain coupling terms were deter- 
mined by fitting to first-principles total energies, vary- 
ing strain and amplitude of uniform LWF distortions in 
the [100], [110] and [111] directions. The interactions 
between LWF coordinates in different unit cells are in- 
cluded to quadratic order only, with the general form 
i ■ l :i'> ■■ Beyond third neighbor, the J ijaf3 

are parameterized as the interaction between two dipoles 
Z*^i and Z*£j, where Z* is the mode effective charge 
for the unstable zone-center phonon, screened by the 
electronic dielectric constant Sac. Z* and are com- 
puted directly using LAPW linear response, while the 
short-range Jij a /3 are fit to the first-principles dynamical- 
matrix elements. 

Using H e ff, classical molecular dynamics simulations 
were carried out for a 10 x 10 x 10 simulation cell, corre- 
sponding to 5000 atoms, with periodic boundary con- 
ditions; the £i a 's in H e ff are in units of 10a, where 
a = 4.016A is the lattice constant. A variable cell 
shape formalism was used together with Nose-Hoover 
thermostats to equilibrate the MD runs at constant tem- 
perature. |H| After equilibration, and prior to computing 
the static and dynamic structure factors, the thermostats 
were turned off and the cell shape and volume were kept 
fixed. Further equilibration (constant-energy MD) gen- 
erally caused the temperature to change by about 5 K. 
After this last equilibration, MD runs of typically 20000 
time steps (each time step ~ 1 femtosecond) were per- 
formed. The static and dynamic structure factors and 
autocorrelation functions of the £j Q 's were then computed 
using data from every 10th time step. 
The different structural phases and transition temper- 
atures T c are identified in the MD simulations by cal- 
culating the three components of the order-parameter 
Sat defined as an average over the LWF coordinates: 
S a = Cia- F° r example, the time average of 

all three components of the order parameter is zero at 
400 K in KNbC>3 , indicating that the system is in the cu- 
bic paraelectric phase. At 370 K, one of the components 
of the order parameter freezes out with a non-zero aver- 
age value of about 0.16, but the time average of the other 
two components is still zero, indicating that the system 
is in the tetragonal phase. Subsequent freezing out of the 



TABLE I. Comparison of calculated and measured transi- 
tion temperatures (see text), between the rhombohedral (R), 
orthorhombic (O), tetragonal (T) and cubic (C) phases. Tem- 
peratures are in Kelvin. 
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a Ref. §). 
b From Ref. [§. 

c See for example, M. D. Fontana, G. Metrat, J. L. Servoin 
and F. Gervais, J. Phys. C 16, 483 (1984). 



other components signals transitions to the orthorhom- 
bic and rhombohedral phases. The MD T c values were 
calculated by us for KNbC>3 and, as a calibration of the 
method, for BaTiOs, using the H e ff constructed in Ref. 
B . The results, given in Table I, are an average of cool- 
ing and heating runs, and we estimate the error in these 
numbers to be about 5 - 10 K. 

The calculated T c for the cubic-tetragonal transition is 
significantly underestimated in both materials, and the 
error is considerably larger in KNbC>3 than in BaTiC>3 
and PbTi0 3 |o|. T c 's for the R-0 and O-T are in bet- 
ter agreement, with the R-0 agreement being the best. 
Possible sources of this quantitative discrepancy include 
the use of the local density approximation, the neglect of 
higher-order coupling to degrees of freedom outside the 
effective Hamiltonian subspace, and the sensitivity of the 
transition to residual inaccuracies in the parametrization 
of strain coupling. In any case, the nontrivial phase se- 
quence and the trend from BaTiOs to KNb03 are cor- 
rectly reproduced, suggesting that the effective Hamilto- 
nian captures the essential behavior of the microscopic 
fluctuations driving the transitions. 

Turning to the soft-mode vibrational frequencies, the 
temperature dependence of the soft-mode dispersion in 
KNb03 is shown in Fig. |l| and compared with the 798 
K cubic phase inelastic neutron scattering measurement 
of Holma and Chen. jl5| The dispersion and tempera- 
ture dependence were determined from the dispersion of 
pronounced peaks in the Fourier transformed autocorre- 
lation function obtained from the MD simulations. The 
experimental data is taken about 100 K above the cubic- 
tetragonal phase transition. The 400 K MD results are 
30 K above the calculated transition temperature. The 
agreement is good below q ~ 0.2, with the theoretical 
TO branch showing greater dispersion for larger q. This 
branch was unusually difficult to measure for q > 0, 
however, the measured peak intensity being extremely 
low and the background unusually high. 
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FIG. 1. Comparison of the soft-mode dispersion along the 
[100] direction. Neutron scattering data from Ref. [15] in the 
cubic phase of KNb03 at 798 K. In the MD simulations, the 
system is in the cubic phase at 400 and 1000 K and in the 
rhombohedral phase at 150 K. 

The most striking feature in Fig. 1 is that the entire 
TO branch softens by about 2 THz on cooling from 1000 
K to 400 K. This behavior differs from the usual soft- 
mode picture of a displacive phase transition, in which 
softening occurs only in the vicinity of the wavevector 
associated with the structural phase transition. It is 
consistent, however, with the previously obtained first- 
principles LAPW linear response results. These re- 
vealed regions of instability in the BZ described by three 
mutually perpendicular interpenetrating slabs centered 
at r, perpendicular to the cubic axes, and extending 
to the BZ boundaries. The first-principles unstable TO 
mode is seen in Ref. |jj to also disperse slightly upwards 
in frequency from T to the X point, consistent with what 
is found here. In the MD simulations, however, the har- 
monically unstable first-principles TO branch is anhar- 
monically stabilized. 

The fact that an entire phonon branch softens, i.e. 
the softening occurs simultaneously over large regions of 
the BZ, suggests that the analysis should focus on the 
temperature dependence of the structure in real space 
rather than in reciprocal space. Tending more to a 
order-disorder picture, this also would naturally incor- 
porate the interpretation of the experimental observa- 
tions of local atomic structure, such as the streak pat- 
terns observed in diffuse X-ray scattering in BaTiOa and 
KNb03. |2j Comes et al. invoked scattering from dis- 
ordered finite-length static chains of distorted primitive 
cells to explain this behavior, proposing an empirical 
model in which there is sequential ordering of chains di- 
rected along the three cubic axes. Thus, for example, 
in the cubic to tetragonal transition, the ordering of the 
z-axis chains corresponds to the disappearence of the in- 
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FIG. 2. Time dependence of the LWF coordinate compo- 
nents and in the orthorhombic (230 K) and cubic (390 
K) phases of KNbC>3, for six (i = 0...5) adjacent primitive cells 
along an uncondensed chain in the ^-direction (see text). R y 
and R z are condensed but R x = 0. Chain-like correlations are 
evident in the correlated motion of the six LWF coordinates 
£i x as a function of time, but not in the y components, which 
oscillate randomly about the condensed value of -0.16. 

coherent scattering (circular streaks in Fig. 6a of Ref. 

due to randomly oriented z-chains. Subsequent or- 
dering of the x-axis and y-axis chains corresponds to en- 
tering the orthorhombic and rhombohedral phases, with 
all the chains being ordered in the ground state rhombo- 
hedral phase. The ab initio linear response calculations 
of Yu and Krakauer [Q on KNb03 provided theoretical 
support for the empirical chain model of Comes et al. by 
showing the existence of BZ planar instabilities. Further- 
more, simulated diffuse elastic X-ray scattering intensi- 
ties calculated from MD simulations on BaTiOs repro- 
duce the observed position and sequential disappearance 
of the three families (circles, vertical and horizontal) of 
streak patterns on cooling from the cubic to rhombohe- 
dral phases. jl6| 

Diffuse X-ray scattering cannot, however, distinguish 
between static and dynamic chains. Using MD simula- 
tions, we can resolve this issue. Fig. || shows the time 
dependence of the LWF coordinates in the orthorhom- 
bic (230 K) and cubic (390 K) phases of KNb0 3 , for 
six adjacent primitive cells along the x-direction. In the 
cubic phase all three order parameters are zero. In the 
orthorhombic phase the order parameters S y and S z are 
condensed, but the average value of S x is zero, corre- 
sponding to disordered uncondensed chains along the x- 
direction. Chain-like correlations are evident in the cor- 
related motion of the longitudinal LWF components as a 
function of time in this figure. (Due to the use of peri- 
odic boundary conditions in the 10 x 10 x 10 simulation 
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FIG. 3. The solid lines are two-site equal-time probabil- 
ity distributions, -P' 2 '(£(i) — £(0)), for i=l,2,3, and 5 along 
condensed y and z chains and for uncondensed x chains in 
orthorhombic KNbC>3. In both cases, £ represents the lon- 
gitudinal component of the LWF coordinate along the chain 
direction. The dashed lines are obtained using the one-site 
probability distribution, assuming no correlation between the 
two sites. The inset shows the one site probability distri- 
butions, P(£(0)), for the x, y, and z components in the or- 
thorhombic phase, as well as in the cubic phase. 

cell, the largest nearest neighbor distance along a chain 
is 5a.) As expected from the linear response calculations 
, the y- and ^-components of the atoms in an x-chain 
should be uncorrelated, and this is confirmed in Fig. |^. 
Similarly if we examine the longitudinal (or transverse) 
LWF components along condensed y- or z-chains (not 
shown) there is no correlation. Chain correlation is still 
evident even as high as 1000 K in the cubic phase, but 
the characteristic chain reversal frequency increases as a 
function of temperature. Fig. || gives a more quantitative 
account of the correlated motions. 

These results unambiguously show 1) the existence in 
real space of chains and 2) the dynamic character of the 
chains. This agrees with the conclusions of Holma et al. 
jl7j , based on their diffuse X-ray measurements, and with 
the empirical lattice dynamics model of Hiiller jl^] . The 
MD simulations show that the chains are preformed well 
above the cubic-tetragonal phase transition temperature. 
The chains are defined by rows of distorted primitive cells 
oriented along the three cubic axes, with the atomic dis- 
placements along the chain highly correlated with one 
other. Displacements in different chains are uncorrelated 
at high temperature, and the observed phase transitions 
correspond to the sequential freezing or onset of coher- 
ence of families of chains along the three cubic axes. The 
softening of entire branches of unstable TO modes gives 
rise to these real space chains and provides a framework 
for understanding both the displacive and order-disorder 
characteristics of these phase transitions. 
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